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Extending the central concept of recurrence times for 
a point process to recurrent events in space-time allows 
us to characterize seismicity as a record breaking process 
using only spatiotemporal relations among events. Link- 
ing record breaking events with edges between nodes in 
a graph generates a complex dynamical network isolated 
from any length, time or magnitude scales set by the 
observer. For Southern California, the network of recur- 
rences reveals new statistical features of seismicity with 
robust scahng laws. The rupture length and its scaling 
with magnitude emerges as a generic measure for distance 
between recurrent events. Further, the relative separa- 
tions for subsequent records in space (or time) form a 
hierarchy with unexpected scaling properties. 



1. Introduction 

Fault systems as the San Andreas fault in California 
or the Sunda megathrust (the great tectonic boundary 
along which the Australian and Indian plates begin their 
descent beneath Southeast Asia) are prime examples of 
self-organizing systems in nature [Rundle et a/., 2002]. 
Such systems are characterized by interacting elements, 
each of which stays quiescent in spite of increasing stress 
acting on it until the stress reaches a trigger threshold 
leading to a rapid discharge or "firing". Since the inter- 
nal state variables evolve in time in response to external 
driving sources and inputs from other elements, the firing 
of an element may in turn trigger a discharge of other el- 
ements. In the context of fault systems, this corresponds 
to earthquakes, or the deformation and sudden rupture 
of parts of the earth's crust driven by convective motion 
in the mantle. 

Fault systems — and driven threshold systems in gen- 
eral — exhibit dynamics that is strongly correlated in 
space and time over many scales. Their complex spa- 
tiotemporal dynamics manifests itself in a number of 
generic, empirical features of earthquake occurrence in- 
cluding clustering, fault traces and epicenter locations 
with fractal statistics, as well as scaling laws like the 
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Omori and Gutenberg-Richter (GR) laws (see e.g. Refs. 
[Turcotte, 1997; Rundle et al, 2003] for a review), giv- 
ing rise to a worldwide debate about their explanation. 
Resolving this dispute could conceivably require measur- 
ing the internal state variables — the stress and strain 
everywhere within the earth along active faults — and 
their exact dynamics. This is (currently) impossible. 
Yet, the associated earthquake patterns are readily ob- 
servable making a statistical approach based on the con- 
cept of spatiotemporal point processes feasible, where the 
description of each earthquake is reduced to its size or 
magnitude, its epicenter and its time of occurrence. De- 
scribing the patterns of seismicity may shed light on the 
fundamental physics since these patterns are emergent 
processes of the underlying many-body nonlinear system. 

Recently, such an approach has brought to light new 
properties of the clustering of seismicity in space and 
time [Bak et a/., 2002; Corral, 2003, 2004; Davidsen and 
Goltz, 2004; Davidsen and Paczuski, 2005; Baiesi and 
Paczuski, 2005], which can potentially be exploited for 
earthquake prediction [Goltz, 2001; Tiampo et aL, 2002; 
Baiesi, 2006]. One aim has been to evaluate distances 
between subsequent events, including temporal and spa- 
tial measures. The observed spatiotemporal clustering of 
seismicity suggests that subsequent events are to a cer- 
tain extent causally related. It further suggests that the 
usual mainshock/aftershock scenario — where each event 
has at most one correlated predecessor — is too simplis- 
tic and that the causal structure of seismicity could ex- 
tend beyond immediately subsequent events, especially 
since the determination of the sequence is largely arbi- 
trary depending on the size of the region considered and 
the completeness of the record of events. 

In this work we quantify the spatiotemporal cluster- 
ing of seismicity in terms of a sparse, directed network, 
where each earthquake is a node in the graph and links 
connect events with their recurrences. This general net- 
work picture allows us to characterize clustering by using 
only the spatiotemporal structure of seismicity, without 
any additional assumptions. 

2. Method 

The key advance we propose is to generalize the no- 
tion of a subsequent event to a record breaking event, one 
which is closer in space than all previous ones, up to that 
time. Consider a pair of events, A and B, occurring at 
times tA < ts- Earthquake 5 is a recurrence of A - or 
record with respect to A - if no intervening earthquake 
hap pens in the spatial disc centered on A with radius 
AB during the time interval [tA,ts]. Each recurrence is 
characterized by the distance / = AB and the time in- 
terval T — ts — tA between the two events. Since the 
spatial window is centered on the first event, any later 
recurrence to it is closer in space than all previous ones, 
and for that reason constitutes another record breaking 
event. ^ This gives rise to a hierarchical cascade of re- 
currences, where each recurrence is, by construction, a 
record. Note that each earthquake has its own sequence 
of records or recurrences that follow it in time. 

Our definition of recurrent events is based solely on 
spatiotemporal relations between events and minimizes 
the influence of the observer by avoiding the use of any 
space, time, or magnitude scales other than those explic- 
itly associated with the earthquake catalog (i.e. its mag- 
nitude, spatial, and temporal ranges). Even the influence 
of the later scales is rather small since, for example, an 



DAVIDSEN, GRASSBERGER AND PACZUSKI: EARTHQUAKE RECURRENCE 



X- 3 



increase in the spatial-temporal coverage of the catalog 
does not generally turn a record-breaking event in a non- 
record breaking event, thus, conserving the property of a 
record. Our definition further allows us to discuss spatial 
and temporal clustering, without introducing any artifi- 
cial scales, or making any arbitrary assumptions about 
the form of seismic correlations. Also, as time goes on, 
one wants to be more strict in declaring B a recurrence of 
A, or related to A in a meaningful way, which is precisely 
what our definition achieves. 

To construct a network we represent each earthquake 
as a node, and each recurrence by a link between pairs 
of nodes, directed according to the time ordering of the 
earthquakes. Distinct events can have different numbers 
of in- going and out-going links, which designate their re- 
lations to the other events. The out-going links from any 
node define the structure of recurrences in its neighbor- 
hood and characterize the spatiotemporal dynamics of 
seismicity, or its clustering with respect to that event. 
The overall structure of the network describes the clus- 
tering of seismic activity in the region that is analyzed. 

To test the suitability and robustness of our method 
to characterize seismicity, we study a "relocated" earth- 
quake catalog from Southern California ^ which has 
improved relative location accuracy within groups of 
similar events, the relative location errors being less 
than 100m [Shearer et al, 2003]. The catalog is as- 
sumed to be homogeneous from January 1984 to De- 
cember 2002 and complete for events larger than mag- 
nitude nic = 2.5 [Wiemer and Wyss, 2000]. Restrict- 
ing ourselves to epicenters located within the rectangle 
{120.5° W, 115.0° W) X (32.5°iV,36.0°iV) and to magni- 
tudes m > rric gives N = 22217 events. In order to 
test for robustness and the dependence on magnitude, 
we analyze this sub-catalog and subsets of it that are ob- 
tained by selecting higher threshold magnitudes, namely 
m = 3.0,3.5,4.0 giving N = 5857,1770,577 events, or 
a shorter period from January 1984 to December 1987 
giving N — 4744 events for m = rric. 

3. Results &6 Discussion 

Fig. 1 shows the probability distribution function 
Pm(l) of distances, /, of recurrent events for different 
thresholds m. The typical or characteristic distance, 
l*{m), where the distribution peaks, increases with mag- 
nitude. For sufficiently large Z, all distributions show a 
power law decay with an exponent ?^ 1.05 up to a cutoff. 
This cutoff is the size of the region of Southern California 
that we consider. 

With a suitable scaling ansatz, the different curves in 
Fig. 1 fall onto a universal curve, except at the cutoff, 
which is a man-made scale imposed on the geological sys- 
tem. The inset in Fig. 1 shows results of a data collapse 
using 

Pm(0-^^"''°'^(V10°''''")- (1) 

The scaling function F has two regimes, a power-law in- 
crease with exponent ~ 2.05 for small arguments and 
a constant regime at large arguments. The transition 
point between the two regimes can be estimated by ex- 
trapolating them and selecting the intersection point, giv- 
ing Lo = 0.012km. For the characteristic distance that 
appears in F we thus find /* ^ Lo x lO^ "^^"^. This is 
close to the estimated behavior of the rupture length 
Lr ^ 0.02 X 10"^/^ km given by Kagan [2002] and re- 
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markably close to Lr = ^/A^ ^ 0.018 x lO^ "^^ ^ km 
given by Wells and Coppersmith [1994], where Ar is the 
rupture area. 

The agreement between our result and that of Wells 
and Coppersmith [1994] suggests that the characteris- 
tic length scale of distances of recurrent events is the 
rupture length, defined in terms of the rupture area 
r = Lr = \/Ar. This is substantially supported by 
the remarkable fact that, for fixed m, Pm{l) and thus /* 
does not significantly vary with the length of the obser- 
vation period despite huge differences in the number of 
earthquakes N — which is very different from a random 
process [Davidsen et a/., 2006]. As Fig. 1 shows, P2.5(0 
is largely unaltered if only the sub-catalog up to 1988 is 
analyzed. This is not true for sub-catalogs of similar size 
generated by randomly deleting events. The comparison 
of the two different observation periods in Fig. 1 further 
shows that /* does not depend strongly on the total num- 
ber of recurrences (or links) or on the average degree of 
the network, (k) = #links/#nodes {{k) = 6.56 (7.40) for 
events up to 1988 (2002) and m — 2.5), but clearly on m. 
The independence of the time span and consequent num- 
ber of events implies that Eq. (1) is a robust, empirical 
result for seismicity. 

The identification Z* = Lr is also consistent with the 
fact that the description of earthquakes as a point pro- 
cess breaks down at the rupture length. Below that scale, 
the relevant distance (s) between earthquakes is not given 
solely by their epicenters but also by the relative loca- 
tion and orientation of the spatially extended ruptures. 
Due to different orientations we expect randomness or 
lack of correlations between epicenters for distances be- 
low the rupture length. If events are happening randomly 
in space, or are recorded as happening randomly in space 
due to location errors, then Pm{l) rises linearly. To see 
this consider a two dimensional disc of radius i?, with one 
point at the center and Nr randomly distributed points. 
The probability that there will be no (other) point within 
a distance / of the center point is (1 — /R^)^^] there- 
fore, the probability density for the closest point to be 
at distance I is {2NrI/R^){1 - l^/R^)^^-^. At smafi Z, 
this will describe the distribution shown in Fig. 1 and 
determine the scaling function F in Eq. (1). In fact, 
this is precisely what the earthquake data show for dis- 
tances smaller than the rupture length (see the straight 
line with a slope of 2.05 in the inset of Fig. 1 and the 
linear increase with slope 1 in the main part of Fig. 1). 

The lengths /* observed for the values of m we con- 
sider are larger than the length 100m) at which we 
observe random behavior due to location errors. In fact, 
the data do not show any anomaly near 100m. More- 
over, Pa{1) (blue triangles) does not change substantially 
if the epicenters in the catalog are randomly relocated 
by a small distance up to one kilometer. Yet, the maxi- 
mum for P2.5(0 shifts to larger / with this procedure, de- 
stroying the scaling of /*(m). Since the smallest /* that 
obeys the data collapse is ^ 160 m, the data collapse 
we observe for the original data verifies that the relative 
location errors are indeed less than 100m, or of that or- 
der. Furthermore, our observations indicate that spatial 
correlations between epicenters are already lost for dis- 
tances 100m < Z < Z*, although the frequency of pairs 
of recurrent events with these small distances is much 
higher than by random chance [Davidsen et al., 2006]. ^ 

Related to the distribution of distances of recurrent 
events is the distribution of distance ratios h/h-i in the 
cascade of recurrences to a given event. Here recurrences 
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are ordered by time; recurrence i comes after i — 1. We 
take lo = 448.5 km, which is the size of the region cov- 
ered by the catalog (Fig. 2a). By construction these 
ratios are always < 1. We denote by Pi(x) the proba- 
bility density that h/h-i = x for each event that has 
an z*^ recurrence. The data for z = 1 (black circles) 
scale over a wide region as Pi(x) x~^'^ with 5r ~ 0.6 
- as already shown in [Davidsen and Paczuski, 2005]. 
This is indicated in Fig. 2a by the straight line. Al- 
though each distribution Pi{x) is different, the curves for 
i > 2 also show (more restricted) power law decay com- 
parable to Pi. For h+i/h 1 they also show a peak, 
which becomes more pronounced with increasing i. This 
is due to recurrences occurring at almost the same dis- 
tance. The observed exponent 5r for the power law de- 
cay has a dynamical origin and is not determined by the 
spatial distribution of seismicity [Davidsen and Paczuski, 
2005] : Purely based on the correlation dimension D2 , one 
would expect Pi{x) ^ x^^~^ . For Southern California, 
this gives a growing dependence Pi{x) ^ rather than 
a decaying behavior. Thus, the exponent Sr reflects the 
complex spatiotemporal organization of seismicity. 

A similar analysis can be made for the distribution of 
recurrence times, Pm{T) for different threshold magni- 
tudes m, which is shown in Fig. 3. These distributions 
all decay roughly as with a ^ 0.9 for intermedi- 

ate times as indicated in the inset. The apparent scaling 
region in Fig. 3 shows some curvature, though. Sur- 
prisingly, Pm{T) is independent of m and the number 
of events in the considered catalog. This is very different 
from earlier results for waiting time distributions between 
subsequent earthquakes [Bak et a/., 2002; Corral, 2003] 
and reflects a new non-trivial feature of the spatiotempo- 
ral dynamics of seismicity that appears when events other 
than the immediately subsequent ones are considered. 

The relative times between subsequent recurrences in 
the hierarchy can be analyzed in the same way as dis- 
tances were above. Fig. 2 shows distributions of ra- 
tios of the times T^/T^+i for subsequent recurrences to a 
given event. The broadest scaling regime materializes for 
T1/T2, again with exponent 6t ~ 0.6. The distributions 
for larger i follow roughly the same behavior for ratios 
^ 1, but deviate (less strongly than for the spatial data 
in Fig. 2a) when the ratios tend to 1. Again it is ob- 
vious that this behavior cannot be explained by random 
events. 

The description of seismicity as a network of earth- 
quake recurrences allows its characterization by means of 
the usual characteristics that are thought to be impor- 
tant for complex networks [Albert and Barabdsi, 2002]. 
One such network property is its degree distribution. 
Fig. 2b shows the degree distributions for m = 2.5, 
which is compared to a a Poisson distribution with the 
same mean degree {k) = 7.40 (solid line). A Poisson de- 
gree distribution would be expected if earthquakes epi- 
centers were placed randomly in space and time. While 
the in-degree distribution agrees with such a random net- 
work, the out-degree distribution shows significant devi- 
ations. In particular, the network keeps a preponderance 
of nodes with small out-degree as well as an excess of 
nodes with large out-degree compared to a Poisson dis- 
tribution. This effect is independent of magnitude, as 
an analysis of subsets with higher magnitude threshold 
shows. Note, however, that {k) decreases with m, sim- 
ply because the catalog size shrinks with m. In partic- 
ular, we find {k) = 6.24, 5.20, 4.35 for m = 3.0, 3.5, 4.0, 
respectively. The non-trivial behavior of the out-degree 
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distribution implies in particular that the network topol- 
ogy and, thus, the hierarchial cascade of recurrences or 
records captures important information about the spa- 
tiotemporal clustering of seismicity. ^ 

4. Conclusions 

Our analysis shows that the description of seismicity 
by means of recurrences in space-time allows us to char- 
acterize its clustering behavior using only spatiotempo- 
ral relations between events and to identify new, robust 
scaling laws in the pattern of seismic activity. The pairs 
of recurrent events form a complex network with non- 
trivial statistics. The method allows us to detect the 
rupture length and its scaling with magnitude directly 
from earthquake catalogs without making any assump- 
tions. Our results for the distributions of relative sepa- 
rations for the next recurrence in space and time should 
also have implications for seismic hazard assessment. Fi- 
nally, our findings provide detailed, benchmark tests for 
models of seismicity. 
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Notes 

1. Notice the difference to the definition of an e-recurrence, 
where any event B is considered a recurrence of A if it 
occurs at a spatial distance less than some fixed threshold e 
[Eckmann et al, 1987]. In our definition we do not impose 
any threshold but allow the sequence of events themselves 
to determine which events are recurrences to other ones. 

2. http:/ /www. data.scec.org/ftp/catalogs/SHLK/ 

3. Note that a systematic dependence of the location error on 
magnitude has not been reported in the literature and is 
also not present in the catalog at hand. It is unlikely that 
the characteristic length we see (/*) is merely an artifact 
due to location error growing with magnitude. 

4. Our results are robust with respect to modifications of 
the rules used to construct the network, e.g., using spatial 
neighborhoods such that the construction becomes symmet- 
ric under time reversal or taking into account magnitudes. 
All such modifications have the drawback that they do not 
define a record breaking process consisting of recurrences to 
each event. Our results are also unaltered if we exclude links 
with propagation velocities larger than 6km /sec 0.1% of 
all links). 
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Figure 1. Distribution of distances / of recurrent events 
for sets with different magnitude thresholds m. The dis- 
tribution for m = 2.5 up to 1988 is also shown. Filled 
symbols correspond to distances below 100 m and are 
unreliable due to location errors. The inset shows a 
data collapse, obtained by rescaling distances and dis- 
tributions according to Eq. 1. The full straight line 
has slope 2.05; the vertical dashed line indicates the pre- 
factor Lo in the scaling law for the characteristic distance, 
/* =Lo X lO^-^^^. 
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Figure 2. (a) Distribution of recurrence distance ratios U+i/U. The straight line corresponds to a decay with exponent 
0.6. (b) Distributions of in- and out-degrees of the network for m — 2.5. The given error bars correspond to y^lv{k). (c) 
Distribution of recurrence time ratios Ti/Ti^i. The straight line has slope -0.62. 
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Figure 3. Distributions of recurrence times for different 
threshold magnitudes m. The distribution for m = 2.5 
up to 1988 is also shown. Filled symbols correspond to 
times below 90 seconds which are underestimated and 
unreliable due to measurement restrictions. The inset 
shows the rescaled distributions. 



